Towards finite density QCD with Taylor expansions 
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Convergence properties of Taylor expansions of observables, which are also used in lattice QCD 
calculations at non-zero chemical potential, are analyzed in an effective Nf —2 + 1 flavor Polyakov- 
quark-meson model. A recently developed algorithmic technique allows the calculation of higher- 
order Taylor expansion coefficients in functional approaches. This novel technique is for the first 
time applied to an effective Nf =2 + 1 flavor Polyakov-quark-meson model and the findings are 
compared with the full model solution at finite densities. The results are used to discuss prospects 
for locating the QCD phase boundary and a possible critical endpoint in the phase diagram. 
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I. INTRODUCTION 

Exploring chiral and dcconfining properties of 
strongly-interacting matter at high temperature is one 
of the most interesting issues in modern physics (for re- 
views, see e.g. [1-3]). In particular, the possible ex- 
istence of a critical endpoint (CEP) in the QCD phase 
diagram at non-zero net baryon number density and its 
consequences for the phase structure of QCD at lower 
temperatures are much under debate, e.g. [4]. High en- 
ergy heavy-ion collision experiments at RHIC and SPS 
started to look for experimental evidence of the CEP 
and experiments at new heavy-ion facilities (FAIR and 
NICA), have been designed to probe the relevant high 
density region in the QCD phase diagram. 

Gaining insight into properties of QCD at non-zero, 
even small values of the quark chemical potential (/i) is 
thus of great interest. However, in theoretical calcula- 
tions the relevant region in the QCD phase diagram is 
not easy to access directly in QCD. Model calculations 
or non-perturbative lattice QCD calculations, based on 
approximation schemes, are currently used. 

At vanishing chemical potential the lattice QCD cal- 
culations get refined systematically and become possible 
with almost realistic quark mass values. They indicate a 
temperature-driven crossover transition from hadrons to 
the quark-gluon-plasma phase, for a recent review see [5]. 

However, even in this case the quantitative determina- 
tion of the transition temperature is still a difficult task. 
A careful analysis of cut-off effects is needed to extract 
this quantity in the continuum limit [6, 7]. 
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At non-zero values of the chemical potential the noto- 
rious fermion sign problem prohibits straightforward lat- 
tice QCD simulations. In order to enter at least the small 
chemical potential region of the QCD phase diagram sev- 
eral extrapolation methods such as the rewcighting tech- 
nique, analytic continuation to imaginary chemical po- 
tential or Taylor series expansions have been proposed. 
All these extrapolation methods have their own, intrinsic 
problems and their reliability is still under investigation, 
see e.g. [8-10] for reviews. Under these circumstances, 
it may be useful to apply approximation schemes, used 
currently in lattice QCD calculations, also in model cal- 
culations. Through a comparison with analytic results 
obtainable in QCD-like models, insight may be gained 
for the application of approximation schemes in lattice 
QCD calculations. 

In this work we focus on the Taylor expansion tech- 
nique where a thermodynamic quantity such as the pres- 
sure is expanded in powers of the chemical potential 
\x around vanishing chemical potential, i.e., around the 
point at which Monte Carlo simulations are not hindered 
by the sign problem. This method was developed for 
studies of QCD thermodynamics with the aim to gain 
information also on the location of a possible critical 
point in the QCD phase diagram [11, 12]. In general, the 
relevant Taylor coefficients can be calculated with stan- 
dard simulation techniques used also at /i = [13-15]. 
With this approach the calculation of several thermo- 
dynamic quantities has been extended to non-zero val- 
ues of the chemical potential; results have been com- 
pared with hadron resonance gas model calculations at 
low temperature [16] and quark-gluon gas thermodynam- 
ics at high temperature. Attempts have also been made 
to determine the phase boundary and to infer a possi- 
ble critical endpoint in the phase diagram from conver- 
gence properties of the Taylor series for thermodynamic 
quantities, e.g., the pressure or quark number suscepti- 
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bility [9, 11, 12, 17]. 

So far, basic features of the Taylor series such as its 
convergence properties have received only little atten- 
tion because a proper investigation of these aspects re- 
quires the knowledge of higher-order Taylor coefficients. 
In most QCD studies these coefficients are difficult to ob- 
tain and a systematic analysis of convergence properties 
is not yet possible. At present any attempts to identify 
some precursor effects of a possible critical endpoint thus 
have to rely on information extracted from a few low or- 
der expansion coefficients. 

Recently, a new algorithmic differentiation technique 
has been developed that allows for the calculation of 
higher-order Taylor coefficients in functional approaches 
[18]. This enables one to investigate general convergence 
properties of the Taylor expansion method at least in 
model calculations. Their consequences and prospects 
for the determination of the phase boundary and a pos- 
sible critical endpoint are discussed here in the framework 
of the effective Nf = 2 + 1 flavor Polyakov-quark-meson 
(PQM) model. The model is studied in mean-field ap- 
proximation and the findings from the Taylor expansion 
are compared with the corresponding model mean-field 
solution at finite density [19]. 

The outline of this work is as follows: In Sec. II the 
Taylor expansion method and its application to the anal- 
ysis of different approximants for radii of convergence are 
introduced. We also discuss the improvement of Taylor 
series expansions by a Pade resummation. Basic strate- 
gies for the determination of the phase boundary and the 
location of a critical endpoint in the phase diagram from 
Taylor expansion coefficients are discussed in Sec. III. 
Sec. IV is devoted to the presentation of the model anal- 
ysis where higher-order Taylor coefficients are used to 
locate the phase boundary and the critical endpoint at 
finite chemical potential. We conclude in Sec. V with a 
brief summary and discussion of our findings. 



around fi = 0: 



n=0 



(1) 



The Taylor expansion coefficients are determined by 
derivatives with respect to the chemical potential eval- 
uated at vanishing fx, 
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Due to the CP-symmetry of QCD, the partition function 
is symmetric in /i, i.e. Z{p) = Z(—fi) which reflects 
the invariance of 2 under exchange of particles and anti- 
particles. As a consequence, all odd expansion coeffi- 
cients vanish and the series is even in p/T '. In QCD, 
various quark chemical potentials have to be taken into 
account in this expansion. However, in order to simplify 
our discussion in the following, we focus on one uniform 
quark chemical potential p, = p q = (J-b/S which will not 
restrict our conclusions. 

Besides the pressure, any thermodynamic observables 
like, e.g., the energy density (e), the trace anomaly 
(e — 3p) or the quark number susceptibility Xq can be 
expanded in a similar series with appropriately modi- 
fied Taylor expansion coefficients [17, 21]. Of course, all 
these modified Taylor expansion coefficients are related 
to the initial expansion coefficients of the pressure. An 
instructive example is the Taylor expansion of the quark 
number susceptibility, Xqi which is defined as the second 
derivative of the pressure w.r.t. to fi 



(3) 



This yields the following relation to the corresponding 
Taylor expansion coefficients c n (T) of the pressure 



II. EXTRAPOLATION TO FINITE /j 



J>2 



= gn(n-l)c B (T) (£ 



(4) 



2 (T) 



The Taylor expansion method has been used as an ap- 
proximation scheme to gain information on QCD ther- 
modynamics at high temperature and for small, non- 
vanishing values of the chemical potential. Here, we start 
with a description of some generic model-independent 
properties of the Taylor expansion of, e.g., the QCD par- 
tition function (see also [20]). 



A. Taylor expansion 

At fixed temperature the logarithm of the QCD par- 
tition function, i.e. the pressure p = {T /V)\a.Z, 
may be expanded in a Taylor series in powers of \xjT 



The expansion of other thermodynamic quantities like, 
for instance, the trace anomaly requires additional 
derivatives of the Taylor coefficients w.r.t. the tempera- 
ture. 

In the following we compare the Taylor expansion of 
the pressure and the quark number susceptibility at non- 
zero chemical potential. If a critical endpoint exists in the 
QCD phase diagram then it is expected that the quark 
number susceptibility diverges at that point while the 
pressure stays finite and is still continuous. Although 
the Taylor series for both quantities have identical radii 
of convergence they thus will behave differently on the 
phase boundary. A comparison of both series as the trun- 
cation order is increased might therefore be helpful to 
gain new insights in locating a possible CEP in the QCD 
phase diagram. 
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B. Radius of convergence 

Whether or not the radius of convergence of the Taylor 
expansion of the thermodynamic potential around fi = 
is related to a phase transition in the physical system 
crucially depends on the location of the singularity in the 
complex /i-plane that causes the breakdown of the series 
expansion. Only if this singularity lies on the real axis 
the radius of convergence corresponds to the location of 
a critical point. Information on the location of the singu- 
larity in the complex ^-plane can be deduced from prop- 
erties of the Taylor expansion coefficients themselves. In 
any case, the singularity closest to fj, = will determine 
the radius of convergence of the Taylor series. This can 
be obtained from the behavior of the Taylor expansion 
coefficients of the pressure, 



r = lim 7'2„ = lim 

n— >OC 71— >oo 



C'2n 



C2n+2 



1/2 



(5) 



Other definitions of the radius of convergence are known 
and their application in the context of QCD have been 
discussed in the literature, see e.g. [12]. 

Of course, in almost all applications only a finite num- 
ber of expansion coefficients are known and the limit 
n — > oo cannot be performed. In some cases it might 
be possible to extrapolate the T2n to infinity, making use 
of their expected asymptotic behavior. Otherwise, one 
will have to assume that the radius of convergence r has 
been reached when subsequent estimators, r2„, do not 
change within errors anymore. 

At a given order n the Taylor expansions of different 
observables will give different estimators for the radius 
of convergence, which are all related. For example, the 
estimator for the radius of convergence obtained from 
the quark number susceptibility, r; 
one obtained from the pressure series via 



2„, is related to the 



' 2n 



c 2n 


1/2 ^ 


c 2n+2 





(2n + 2)(2n+ 1) 
(2n + 3)(2n + 4) 



1/2 



f2n+2 ■ (6) 



In the limit n — > oo both estimators for the radius of 
convergence will converge to the same unique limit r since 
the prefactor in front of rm+2 tends to one. However, at 
finite order n the estimates for rm and r\ n will differ 
and as a consequence both may approach the limit r at 
a different rate. For large n we find from Eq. (6) 



"2n = r ?n+2 



1 



1 



+ 0(n~ 2 ) 



(7) 



When deviations from the asymptotic value r are of 
C(l/n), i.e., r2„+2 — r(l + A/n) with a positive, con- 
stant A, both estimators for the radius of convergence 
will converge with the same truncation errors of 0(l/n) 
(cf. App. A). Deviations from the asymptotic value r 
will, however, be numerically smaller in the susceptibility 
series than in the pressure series by a factor (A—l)/A. In 
other words, the estimator for the radius of convergence 



obtained from the susceptibility series at given order n 
might be more suitable for estimating the true radius of 
convergence than the one obtained from other observ- 
ables, e.g., from the pressure series. 



C. Resummations: Pade approximation 

The convergence of a Taylor scries can be improved 
further by a rcsummation of the expansion coefficients 
which is based on a Padc approximation (for details sec 
[22]). For clarity and completeness, we collect the needed 
facts here. The Pade resummation employs the same 
derivative information as the initial Taylor expansion co- 
efficients but often shows an extended convergence range, 
in particular in the presence of singularities. Further- 
more, a more stable result can be achieved with fewer 
coefficients. 

The Padc approximant can be obtained by rewriting a 
Taylor series of an analytic function t(x) around xq = 
to order N 



N 



(8) 



as a ratio of two power series p(x) and q(x) of order L 
and M, respectively 



p(x) _ Po + Pix j h p L x L 

q(x) 1 + q\x H h qnx M 



[L/M] = R LM (x) 

(9) 

One advantage of this reformulation is that poles and 
singularities are better characterized by a rational func- 
tion than by a polynomial Taylor series. The expan- 
sion coefficients for p and q can be obtained up to order 
N = L + M by equating the derivatives of both series up 
to order N at x = 0, i.e., 



3 l i?L,M(z) 



dx l 



dH(x) 



x=0 



dx l 



for i < N . 



(10) 



x=0 



An alternative and more elegant form relates the Padc 
approximant to determinants of two (M + 1) x (M + 1) 
matrices as follows 



[L/M] 





M+l) 


t(L-M+2) 


■ ■■ t 
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t(L-M+3) 
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t(L) 
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t(L+l) 
-M+l 

E %)* MH 

j=0 


■ ■ ■ t{L+M) 
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t(L-M+2) 


t(L-M+3) 


■ ■ t (L+2) 






t(L) 

x M 


t(L+l) 
x"- 1 


' ' t (L+M) 
1 





(11) 

Note that for M = and L = N the Pade approxi- 
mant is identical to the initial Taylor approximation, i.e., 
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[N/0] = t(x). Of course, apart from choosing the highest 
derivative order N in the Pade approximation one can 
also vary L or M independently. 

The determination of the radius of convergence or 
the phase boundary with the Pade approximation is 
more involved. For the case [N/2] one immediately 
sees, Eq. (11), that the Pade approximant has a pole at 
x = ±y/ 'cpf I 'cjv+2 since the odd Taylor coefficients van- 
ish. This pole coincides with the estimator for the radius 
of convergence, r/v, of the Taylor series up to order x N+2 , 
Eq. (5). For a general Pade approximant [L/M] the pole 
structure is more elaborated. However, in this case the 
first pole at positive and real chemical potential might 
provide a useful estimate for the phase boundary. 



III. LOCATING A CRITICAL ENDPOINT 

It has been argued that the estimators for the radius 
of convergence can be used to determine the location of a 
possible CEP in the QCD phase diagram [11, 12, 23]. In 
general, a breakdown of a Taylor expansion is caused by 
a singularity in the complex /i-plane. Only if this singu- 
larity, which is closest to the origin at /j, = 0, lies on the 
real /x-axis, it is also related to a physical phase transi- 
tion since it is a singularity in a thermodynamic quantity 
corresponding to a zero in the partition function [24]. 
If a CEP exists in the QCD phase diagram, it belongs 
to the universality class of the three dimensional Ising 
model. The corresponding second-order phase transition 
will lead to divergences in thermodynamic observables 
with a power law behavior. 

The finite radius of convergence arises from a singu- 
larity on the real /x-axis only if there exists a no such 
that for all n > no all expansion coefficients are positive. 
Since in practice only a few terms in the Taylor series are 
known such a mathematically rigorous conclusion cannot 
be drawn in general. The hope is that by investigating 
the structure of a few known coefficients some hints on 
the properties of the Taylor series and an estimate for 
the location of a CEP in the phase diagram can be ob- 
tained. Through an inspection of the relative magnitude 
and relative signs of the expansion coefficients the deter- 
mination of the radius of convergence of the Taylor series 
might be possible and the singularity in the complex u- 
plane which causes the breakdown of this expansion can 
be located. 

Since the thermodynamics in the vicinity of a second- 
order phase transition is dominated by the singular con- 
tribution of the free energy density or thermodynamic 
potential a generic structure of the Taylor expansion co- 
efficients should occur [25]. At vanishing chemical po- 
tential the singular behavior near a phase transition is 
controlled by the reduced temperature which depends 
linearly on T — T c but quadratically on the chemical 
potential. Thus, a (2n) th derivative with respect to fj, 
shows a singular behavior that is similar to that of the 
n th derivative with respect to temperature. This infers 



that the structure of the coefficient C2 n +2 can be esti- 
mated by investigating the temperature derivative of C2„. 
This is indeed found in our model calculations. For an 
example see Fig. 1 and the discussion in Sec. IV B. This 
means that C2 n +2 w iU stay positive only up to a tem- 
perature which corresponds to the first maximum of the 
previous coefficient C2 n - Near the transition tempera- 
ture the coefficients start to oscillate. Thus, through a 
determination of the first maximum in C2«, an estimate 
for the largest temperature below which all coefficients 
are positive 1 , becomes feasible. Only in this case where 
all coefficients are positive the corresponding singularity 
lies on the real axis. Furthermore, this gives an esti- 
mate for the critical temperature where a possible CEP 
in the phase diagram might be located. Once this critical 
temperature has been obtained the corresponding critical 
chemical potential follows immediately via \i c = r(T c ). 

At higher temperatures the sign of the Taylor coeffi- 
cients has no generic structure and the singularity that 
limits the radius of convergence of the Taylor series will 
move in the complex /x-plane. In this case the conver- 
gence radius provides an upper bound since the crossover 
is related to the real part of the singularity. 

For temperatures below the CEP and larger chemical 
potentials a first-order phase transition emerges. As a 
consequence, two degenerate minima of the grand po- 
tential exist at criticality. Since the Taylor expansion 
method involves only local information in the vicinity of 
one minimum of the effective potential, it is, in general, 
not possible to resolve a first-order phase transition with 
this technique. However, it still should be possible to give 
an upper bound for the location of the phase boundary. 



IV. A MODEL ANALYSIS 

In this section we apply the Taylor expansion method 
to an effective Nf = 2 + 1 flavor Polyakov-quark-meson 
model. This model exhibits a critical endpoint in the 
(T, fi) phase diagram and enables us to inspect the con- 
vergence properties of the Taylor series for first-, second- 
order and crossover phase boundaries. For the first time, 
the Taylor expansion coefficients are calculated to very 
high orders, so that the convergence of the series can be 
probed reliably. Furthermore, the results are confronted 
with a full model solution. 



1 As pointed out earlier the mathematically more rigorous criterion 
for a singularity on the real axis does not rule out the possibility 
to have a few negative expansion coefficients. We assume here 
that once a negative coefficient occurs this will lead to irregular 
signs of the coefficients also in higher orders. 
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A. The Polyakov-quark-meson model 

The three flavor Polyakov-quark-meson model [26] is 
a straightforward extension of the two flavor model [27] 
and serves as an effective model for strongly-interacting 
matter which incorporates a chiral phase transition, sig- 
naled, for example, by a peak in the chiral susceptibility. 
Through the coupling to the Polyakov loop, which is an 
order parameter for deconfinement in the infinite quark 
mass limit, it is also sensitive to deconfining aspects of 
the QCD phase transition. The thermodynamics of the 
PQM model is, for vanishing chemical potential evalu- 
ated in a mean-field approximation, in good agreement 
with recent 2 + 1 flavor lattice data [19, 26, 28]. In con- 
trast to lattice simulations, the model can be applied 
directly at finite chemical potential. 

The PQM model Lagrangian reads 

Cpqm =q{Hp-h<t>$) q + £ m -«(*,$) , (12) 

where q = (u, d, s) denotes the three quark fields. The 
interaction between the quarks and the mesons is imple- 
mented by a flavor-blind Yukawa coupling h. The meson 
matrix cf>5 = X)a=o(^ a /^) ("a + ilfo n a) consists of nine 
scalar a a and nine pseudo-scalar n a meson fields. A a are 
the usual Gell-Mann matrices. The covariant derivative 
Ip= $ — i^foAo couples the Polyakov-loop and its conju- 
gate to the fermionic degrees of freedom where a gauge 
coupling has been absorbed in the gauge fields. For de- 



tails see [27]. 

The purely mesonic contribution is given by 

Cm = Tr (d^8^) - m 2 Tr(^0) - A x [Tr{4f </>)] * 
-A 2 Tr (0V) 2 + c (dct(0) + det(</) t )) 
+ Tr[H(4> + ^)] , (13) 

with 4> = ^ n (A„/2) (a a + iir a ). Chiral symmetry is ex- 
plicitly broken by the last term in Eq. (13) and the 
£/(l)A-symmetry by the 't Hooft determinant term with 
a constant strength c. In the 2 + 1 flavor scenario only 
two order parameters (ao) and (as) emerge in the singlet- 
octet basis. They are conveniently rotated to the non- 
strange (a x ) and strange (cry) basis, see [29] for details. 

The grand potential evaluated in mean-field approxi- 
mation is a sum over the mesonic U, fermionic fl qq and 
Polyakov loop U contributions: 

n = U {(a x ) , (a y )) + n qq ((a x ) , (a y ) *) , 

(14) 

where <& and $ denotes the real Polyakov loop expecta- 
tion values. Note, that at finite chemical potential the 
$ and $ are not linked by complex conjugation and the 
effective Polyakov loop potential depends on two inde- 
pendent variables. 

Explicitly, the quark-anti-quark contribution in pres- 
ence of the Polyakov loop reads 



*.»„>,*,*) =-2T £ /(0a{ ln 



f—u,d,s ' 



In 



1 + 3($ + |. e -( i ^/-A<)/T) e -(£, 3 , / -M/T + e -3(E qif - M )/T 

1 + 3(# + ^ e -^f + ^y T )e-^ E o-' +tl / T + e -3(E qJ+f ,)/T | ( 15 ) 



with the flavor-dependent single-particle energies 



E qJ = Jk 



(16) 



and quark masses 



mi = h (a x ) /2 and m s = h (a y ) /v2 (17) 

for the light (/ = u, d) and strange quarks, respectively. 
In a mean-field approximation, fluctuations are neglected 
and a divergent vacuum contribution to the grand po- 
tential has been omitted here which is irrelevant for our 
following discussion. See [30] for an investigation of the 
effect of this vacuum term. 

The Yukawa coupling h is fixed such to reproduce a 
light constituent quark mass of order m; ps 300 MeV. 
The parameters of the mesonic potential U are fitted to 
well-known pseudo-scalar meson masses and decay con- 
stants [29] . As in [26] we will employ m a = 600 MeV. 



For the Polyakov-loop part several potentials can be 
found in the literature which describe the thermodynam- 
ics equally well in the vicinity of the chiral transition at 
H = 0, see e.g. [26]. Here we employ the logarithmic 
version [31] 



a(T) 
2 

-b(T) In 



(18) 



1 - 6$$ + 4 ($ 3 + $ 3 ) - 3 ($$)' 



with the temperature-dependent coefficients 

2 



( y 



and the parameter To = 270 MeV. Note that we ignore 
any Nf- and/or //-modifications of this parameter here, 
i.e., the matter back-reaction to the gluon sector, in order 
to simplify our discussion, sec [26, 27, 32] for more details. 
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Finally, the temperature and /i-dcpendent order pa- 
rameters are determined by the solution of the gap equa- 
tions which minimize the grand potential 



da x 



on 



dfl _ 90 
5$ _ d¥ 







(19) 



where min labels the (T, /^-dependent global minimum 
of the potential. 



B. Higher-orders Taylor coefficients 



1-10° 



-1-10° 



r 




10c 10 

101/12! (dc 10 

c 12 


dT)7tT 



0.99 



1 

T/T v 



1.01 



Since the grand potential or the pressure are known, 
the Taylor expansion coefficients in Eq. (2) can be ob- 
tained by calculating the corresponding derivatives of the 
grand potential evaluated at the physical potential min- 
imum. This tedious task is further hindered by the im- 
plicit /^-dependence of the order parameters, introduced 
by Eq. (19). The explicit /^-dependence is solely given 
by the fermionic part of the potential Clq q whereof the 
/i-derivatives can, in principle, be obtained analytically. 
However, a fully analytical approach is infeasible since 
the implicit dependences can only be inverted numeri- 
cally. Standard numerical derivative techniques such as 
divided differences may be used for this task at least for 
the first few derivatives, but definitely fail for higher or- 
ders due to the increasing rounding and truncation errors. 

In order to proceed a novel derivative technique, based 
on algorithmic differentiation (AD), was developed that 
properly allows for the implicit dependencies and further 
avoids numerical errors. The AD technique enables the 
evaluation of derivatives to high orders even in the case of 
only numerically solvable implicit dependences, such as 
in Eq. (19). Furthermore, the derivatives can be obtained 
with extremely high numerical precision, essentially lim- 
ited only by machine precision. Details of this method 
can be found in Ref. [18]. 

Of course, this AD approach is not limited to simple 
Taylor series and can also be applied to more involved 
Taylor expansions of physical quantities which also re- 
quire temperature-derivatives. An example would be the 
Taylor expansion of the trace anomaly (e — 3p). 

The Taylor coefficients c„ for the pressure, Eq. (2) 
have been calculated with this novel AD technique for 
the three flavor PQM model and can be found, up to the 
22 nd order, in [19]. All coefficients for n > 4 show rapid 
oscillations near the transition temperature T x at fi = 0. 
The coefficients become small outside a 5% temperature 
interval around the fi — transition temperature T x , i.e. 
0.95 < T/T x < 1.05. Within this temperature interval 
the amplitudes of the oscillations increase with increas- 
ing order and hence higher orders are important even for 
small expansion parameters, i.e. (p/T) < 1. As an ex- 
ample the coefficient C12 and the properly rescaled tem- 
perature derivative of the coefficient Cio are shown in 
Fig. 1. As already discussed in Sec. Ill it obeys the same 
singular behavior as C12, hence the structure of c\i can 



FIG. 1. The Taylor coefficient C12 and the (rescaled) tem- 
perature derivative of the coefficient cio in the vicinity of the 
critical temperature T x . For comparison the (rescaled) Taylor 
coefficient cio is also shown. 



be estimated by the temperature-derivative of Cio- In- 
terestingly, even the magnitude of the coefficient c\i is 
well estimated by dimensional analysis and compensat- 
ing the different factorials. Some convergence analysis as 
well as detailed extrapolations of the pressure and quark 
number susceptibility Taylor series to finite [i with these 
coefficients can be found in [19]. 



C. Determining the phase boundary 



We start with a calculation of the (T, fi) phase di- 
agram of the PQM model without referring to a Tay- 
lor expansion. At vanishing chemical potential the chi- 
ral crossover, signaled here by a peak in the chiral sus- 
ceptibility occurs at T x ~ 206 McV. At this tempera- 
ture also the Polyakov loop susceptibility has its maxi- 
mum. For non- vanishing chemical potentials the position 
of maxima in both susceptibilities differ. The peak in 
the Polyakov loop susceptibility occurs at higher tem- 
peratures than that in the chiral susceptibility. This 
scenario changes, however, when fluctuations are taken 
into account, e.g., if /i-corrections of the To parameter in 
the Polyakov loop potential are considered, see [27, 32]. 
For small temperatures a first-order chiral phase tran- 
sition is found which terminates at larger temperature 
in a CEP. In mean-field approximation the CEP is lo- 
cated at (T c ,/i c ) ~ (185, 167) MeV which corresponds to 
T c ~ 0.9 T x and /i c /T c ~ 0.9. Approximations beyond 
mean-field push the location of the CEP off the T-axis. 
For a mean-field and renormalization group comparison 
in a quark-meson model, sec e.g., [33]. In the following 
we compare several methods for the determination of the 
phase boundary with the Taylor coefficients. 
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FIG. 2. Estimates for the radius of convergence obtained from 
T2n (solid lines) and r^ n _ 2 (dashed lines) for different orders 
(2n = 8, 16, 24) of the Taylor expansion. Also shown are the 
phase boundaries for the chiral (red line; dashed: crossover, 
solid: first order) and a deconfinement line identified by the 
peak position of the Polyakov loop susceptibility (yellow line) . 
The black dot indicates the location of the CEP. 
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FIG. 3. Estimate of the phase boundary at T = 190 MeV ~ 
0.92 T x with r2n and r% n and poles in the Pade approximation 
of the pressure and quark number susceptibility as a function 
of the order of the Taylor expansion 2n. The used highest 
coefficient is C2 n +2- The dotted horizontal line indicates the 
phase boundary calculated directly at finite fi. Note also that 
an estimate is possible only for even Pade approximations. 



1. Taylor expansion 



In Fig. 2 wc show the results for two estimators for 
the radius of convergence, obtained from the expansion 
coefficients of the pressure, r-m (solid lines), and of the 
quark number susceptibility, f2„_ 2 (dashed lines), for 
three different truncation orders (2n = 8, 16, 24) in the 
(T, fi) phase diagram. The phase boundary for the chi- 
ral transition and a deconfinement line, corresponding to 
a maximum in the Polyakov loop susceptibility, are also 
displayed. Note that the difference in the radii of conver- 
gence, estimated from the r2„ and r\ n _^ respectively, is 
only caused by the prefactor in Eq. (6) since the same 
Taylor coefficients contribute to both estimators. 

Near T x the sign of the expansion coefficients oscil- 
lates, indicating that the break down of the series expan- 
sion is caused by a singularity in the complex plane away 
from the real axis. The oscillations in the Taylor coef- 
ficients entail oscillations in the radii of convergence at 
small chemical potentials which are hard to see in Fig. 2. 
Hence, for small chemical potentials the radius of conver- 
gence is not suitable to estimate the phase boundary. For 
larger values of the chemical potential and for tempera- 
tures below T x estimates for the radius of convergence 
approach the phase boundary from above with increas- 
ing truncation order n. The observed agreement suggests 
that the singularity limiting the Taylor expansion is close 
to the real /^-axis and the small imaginary part is negligi- 
ble when compared to the error due to the finite number 
of coefficients available for a fast crossover. In particular, 
this is still valid for the region in the vicinity of the CEP. 
We note that for n>10 corrections to the asymptotic 
value r are well described in terms of C(l/n) corrections. 
In fact, with an ansatz r(l + A/n) for the estimators of 
the radius of convergence r2„, the leading C(l/n) cor- 



rections can be eliminate by using subsequent estimates. 
This leads to estimates for r obtained from the pressure 
and quark number susceptibility series that agree within 
a few percent and reach the asymptotic value within 15% 
for 2n > 14. The coefficient A for the correction term, 
obtained in this way, is about a factor 5 larger in the 
pressure series than in the quark number susceptibility 
series. Qualitatively, this different asymptotic behavior 
of both series can be understood in terms of the structure 
of the singular contribution to the grand potential that 
limits the radius of convergence. We discuss this in more 
detail in the Appendix A. 

At low temperatures the estimates for the radius of 
convergence increase with increasing order n of the Tay- 
lor expansion and allow also ratios /i/T > 1. This is 
in agreement with the small magnitude of the higher or- 
der coefficients in this temperature region. In this re- 
gion the results are in agreement with results from the 
resonance gas (cf. [16]), i.e. the estimates of the con- 
vergence radius of the pressure are given by r™ G (T) = 
1/3^(1 + 2n)(2 + 2n). 

As already mentioned above, it is not possible to re- 
solve a first-order transition with the Taylor expansion. 
Therefore, T2 n as well as r\ n yield estimates for the ra- 
dius of convergence that are larger than the true location 
of the first-order transition line and signal a misleading 
convergence of the Taylor expansion. In this region the 
Taylor expansion might provide stable and continuous re- 
sults but one cannot extract the precise position of the 
boundary for a first-order phase transition. 
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2. Pade approximation 

The Pade approximation could improve the conver- 
gence behavior and therefore provide a better estimate of 
the phase boundary since the Pade approximant employs 
all Taylor coefficients in contrast to the previous analysis 
where only two subsequent Taylor coefficients enter. The 
pole in the [N/2] Pade approximation yields the corre- 
sponding radius of convergence as described in Sec. II C. 
For a general [N/M] Pade approximation we use the first 
pole at real and positive /i in order to estimate the phase 
boundary. Since all Taylor coefficients with their corre- 
sponding error enter in the Pade approximant the error 
propagation is more involved here in contrast to the pre- 
vious discussion where only two error sources enter. In 
this context, the application of the Pade approximation 
in lattice simulations is much more involved. However, 
in the model analysis the coefficients, obtained with the 
AD technique, exhibit extremely small numerical errors 
and result in a stable and reliable Pade approximation. 

Apart from the diagonal scheme [N/N] we also use a 
non-diagonal [N + 2/N] approximation for the suscepti- 
bility and a [N + 2/N - 2], respectively [N + 2/N] ap- 
proximation for the pressure and determine always the 
first pole at real and positive fx. 

In Fig. 3 we show for a fixed temperature T = 0.92 T x 
the estimates of the phase boundary extracted from Pade 
approximations for the pressure (p) and the quark num- 
ber susceptibility (x) for different truncation orders n. 
The chosen temperature is slightly above the tempera- 
ture of the CEP in the phase diagram. In addition, the 
estimates for the radii of convergence (r2„ and r^J and 
the chiral crossover chemical potential (dashed horizontal 
line) are also shown. For small truncation orders 2n < 8 
the [A r /2] Pade approximations coincide with the radii 
of convergence as expected. For 2n > 8 we observe that 
the Pade approximation converges faster for the pressure 
as well as for the susceptibility compared to the radii of 
convergence obtained without the Pade approximation. 

We conclude that the Pade approximation converges 
faster at larger truncation orders, in particular for the 
quark number susceptibility. The estimate of the phase 
boundary from the Pade approximation becomes com- 
parable to the one obtained from estimators for the ra- 
dius of convergence determined directly from ratios of 
Taylor expansion coefficients only at significantly larger 
truncation order. For example, in Fig. 3 the distance of 
the Pade estimate to the horizontal line at 2n = 12 is 
achieved with the Taylor expansion only for 2n > 20. 
However, the lower truncation order in the Pade approx- 
imation induces a more involved error propagation. This 
is no problem with the AD-tcchniqucs but might hamper 
lattice simulations. 

a. First-order transition As already mentioned, the 
Taylor expansion around ll = fails if a first-order phase 
transition emerges since the jump to a new global mini- 
mum of the effective potential cannot be resolved. How- 
ever, it might still be possible to enter the metastable 
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FIG. 4. Taylor expansion of the (non-strange) order param- 
eter a x compared to the direct evaluation at finite /i for 
fi/T — 3. The dotted lines indicate the beginning of the 
metastable phase. The vertical line indicates the estimate for 
the radius of convergence, r24- 

phase with the Taylor expansion method. In this case 
the radius of convergence just reflects the fact that this 
expansion is still valid in the metastable phase. However, 
a precise determination of the first-order transition is not 
possible. This conceptual drawback is not modified when 
the Pade approximation is employed. 

As an example we show in Fig. 4 for fi/T = 3 the chi- 
ral order parameter, the non-strange condensate o~ x , as 
a function of the temperature in the vicinity of a first- 
order transition. Here o~ x has been obtained from a Tay- 
lor series truncated at 2A th order. The lines, labeled with 
PQM, represent the full model solution of the gap equa- 
tion Eq. (19) where a first-order transition can be re- 
solved. Both solutions agree very well in the low temper- 
ature phase until T ~ 96 MeV. The dotted line indicates 
the begin of the metastable phase. One sees that it is still 
possible to enter this metastable phase with the Taylor 
expansion method. The vertical line in the figure marks 
the location of the estimate for the radius of convergence 
obtained for the chosen n = 24 truncation order. Fur- 
thermore, this example also demonstrates that the Taylor 
expansion method is applicable also for ll/T > 1 as long 
as one stays within the convergence regime. 

b. Critical endpoint Directly at the critical endpoint 
a second-order phase transition emerge and no disconti- 
nuity in the order parameter appears. Hence, the radius 
of convergence might be well-suited to estimate the loca- 
tion of the critical endpoint in the phase diagram once 
the critical temperature is known. The corresponding 
critical chemical potential can again be obtained via 

Li c = r{T c ) = lim r n {T c ) . (20) 
n—tco 

From the discussion given so far and from Fig. 3 one 
may conclude that at fixed temperature the radius of 
convergence can be estimated with an accuracy of about 
15%-20% from a series of the order C(^t 12 ). This is clearly 
an acceptable uncertainty in view of the difficulties one 
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FIG. 5. Temperature T Cj „ at which the Taylor coefficients c„ 
becomes negative for different truncation orders n. The solid 
horizontal line denotes the critical temperature of the CEP. 

has to face in current QCD calculations with non-zero 
chemical potential. 

However, the determination of T c is still a non-trivial 
task. A criterion is needed, that allows to estimate at 
a given order of the expansion the temperature regime 
in which all expansion coefficients may stay positive as 
discussed in Sec. III. As an example we show the ze- 
ros of the Taylor series of the pressure in the complex 
///T-plane around a temperature point where a coeffi- 
cient changes its sign in Appendix B. 

In Fig. 5 we plot the temperature T c>n defined as the 
temperature where the Taylor coefficient c„ becomes neg- 
ative for different truncation orders. This yields an upper 
bound for the critical temperature. Of course, in the limit 
n — > oo the temperature series of these approximations 
should approach the critical temperature of the CEP, i.e., 

T c = lim T c .„ . (21) 

n— ^oo 

However, at least with our model analysis, only a very 
slow convergence of the estimators for the critical tem- 
perature at the CEP is visible. Even at high truncation 
order of 20 th there is still a temperature difference of 
the order of 20 MeV which is why we do not reconsider 
Eq. (20). 

In addition, this convergence behavior might be im- 
proved by applying the Pade approximation but the iden- 
tification of a proper convergence criterion is more in- 
volved as discussed in Appendix B. 

V. SUMMARY AND CONCLUSIONS 

The convergence properties of the Taylor expansion 
method which is used to extend lattice QCD calculations 
to non-zero chemical potential have been investigated. 
The method has been applied to an effective Nf = 2 + 1 
flavor Polyakov-quark-meson model and the results ob- 
tained with the Taylor expansion have been compared to 



a full model mean-field solution. By means of a novel 
algorithmic differentiation technique higher-order Taylor 
coefficients could be calculated for the first time. As a 
function of the temperature the coefficients start to os- 
cillate in the vicinity of the phase transition region. The 
number of roots and magnitude of the oscillations in- 
crease with the coefficient order. Two definitions of the 
radius of convergence, one through the pressure series 
coefficients and the other one via the corresponding sus- 
ceptibility coefficients have been discussed. The conver- 
gence of the estimators for the radius of convergence of 
the susceptibility series is faster than the corresponding 
one for the pressure series. 

With the Pade approximation the convergence rate 
could be further improved. By investigating the pole 
structure of the Pade approximant a definition for the 
radius of convergence could be extracted and a better 
estimate for the phase boundary could be achieved. 

While the convergence radius yields the precise loca- 
tion of the phase boundary only for a second-order tran- 
sition, i.e., at the critical endpoint, we observed a good 
agreement of the convergence radius with the full model 
solution also for temperatures above T c . In view of the 
error of convergence radius due to the finite number of 
coefficients the convergence radius can provide a valuable 
estimate also for the location of a rapid crossover where 
the limiting singularity is close to the real axis in the 
complex /i-planc. 

The singularity lies on the real axis only in the case 
where all coefficients are positive. This criterion also 
provides an estimate for the critical temperature T c . n . 
In our study this estimate depends only weakly on the 
truncation order. But even for truncation orders above 
20 th there is still a significant gap to the full model so- 
lution and a precise determination of the location of the 
CEP remains open. Note that the error in estimating T c 
also affects the estimation of \x c . 

Although, the situation may be different in model cal- 
culations and in QCD, our study suggests that for a 
meaningful application of the Taylor method to lattice 
QCD higher order Taylor coefficients are definitely nec- 
essary in order to control systematically truncation er- 
rors in the series expansion and to estimate reliably the 
location of a possible critical endpoint in the QCD phase 
diagram. At present the available number of coefficients 
extracted in lattice QCD calculations is not sufficient to 
locate a critical endpoint in the QCD phase diagram with 
the Taylor expansion method. However, if the critical 
endpoint is located at smaller temperatures a broader 
temperature interval is expected to emerge in which the 
Taylor coefficients the Taylor coefficients will oscillate. 
The estimators for the critical temperatures T Cj „ are then 
expected to show a stronger temperature dependence and 
they should decrease faster. However, in this case also the 
CEP may be located at large values of [ijT and higher 
orders in the expansion may be needed due to this. 

Furthermore, when going beyond mean-field approxi- 
mations in our model calculation, quantum and thermal 
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fluctuations in the hadronic phase will surely modify the 
oscillations of the Taylor coefficients. This too may im- 
prove the estimate of the critical temperature with the 
Taylor expansion method. However, also in this case the 
location of a possible CEP is shifted towards the /x-axis 
[34] and more coefficients are necessary for an adequate 
description. To explore this quantitatively we plan to 
extend the present study beyond mean-field as well as 
tightening the AD techniques with functional rcnormal- 
ization group methods. 

With the Taylor expansion method it is in general not 
possible to describe a first-order transition completely in- 
cluding the precise determination of the critical tempera- 
ture since only information about one potential minimum 
are available. One may, however, establish upper bounds 
for the location of the phase boundary. 
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Appendix A: Asymptotic behavior of Taylor 
expansion coefficients 

The QCD critical point is expected to belong to the 
universality class of the three dimensional Ising model. 
The critical behavior in the vicinity of this point is con- 
trolled by two couplings which characterize fluctuations 
in thermal and magnetic (symmetry breaking) directions 
of the effective Ising Hamiltonian. The chemical poten- 
tial will mix with both couplings [35] . When approaching 
the critical point in the (T, /j,)-plane the dominant singu- 
lar contributions to, e.g., the quark number susceptibil- 
ity arise from the magnetic coupling. When approaching 
the critical endpoint in the QCD phase diagram at fixed 
temperature, the quark number susceptibility thus is ex- 
pected to scale like [33, 35, 36] 

Xq (T, n)/T 2 ~ (1 - m/Mc) 17 ^ 1 = (1 - v/^V /f)5 (Al) 

which arises from a singular term in the grand potential, 
i.e., the pressure which has the form, 

(p/T 4 ) sing ~(l-/x/M c ) 1+ ^. (A2) 

Assuming that these singular terms give indeed the dom- 
inant contribution to the Taylor expansion coefficients at 



high order, we obtain for the asymptotic behavior of the 
estimators for the radius of convergence (r = /i c ) of the 
pressure and quark number susceptibility series, respec- 
tively, 

/ 25+l\ 

^~Mc(l + ^). (A3) 

We thus expect that in comparison to the pressure se- 
ries, corrections to the radius of convergence are smaller 
in the quark number susceptibility series. In mean-field 
approximation, which is used in our current study, the 
relative factor is 28 + 1 = 7 and it is about 10.6 in the 
three dimensional Ising model. Indeed, this seems to be 
realized with the series discussed in this study, although 
the approach to the asymptotic value is slow. 



Appendix B: Roots in the complex /i/T-plane 

As mentioned in Sec. Ill a singularity in the complex 
/x/T-plane limits the radius of convergence of a Taylor 
expansion. In Fig. 6 an example of the root distribution 
in the complex /.t/T-plane of a 16 th order Taylor series 
of the pressure is given. The temperatures in the three 
panels are chosen such that the highest Taylor coefficient 
Ci6 changes sign while all the remaining coefficients stay 
positive. In the left panel c\q is still positive. As ex- 
pected all 16 roots are complex and distributed on a cir- 
cle whose radius corresponds to the radius of convergence 
ri6- For Ci6 = (middle panel) the Taylor polynomial 
is of degree 14 and consequently only 14 roots emerge. 
Approaching ciq = from positive values the two root 
pairs which are each closest to the imaginary axis merge 
to each one purely imaginary root and the other root van- 
ishes at imaginary plus or minus infinity. When the tem- 
perature is further increased, i.e., c\q becomes negative 
(right panel), two roots reenter the complex zx/T-planc 
on the real /i-axis from plus or minus infinity and rapidly 
approach the radius of convergence. 

A similar analysis could also be performed with the 
Pade approximation which shows, in general, a more 
rapid convergence behavior. However, the investigation 
for this case is more involved since two polynomials oc- 
cur and possible root cancellations can emerge which are 
hard to distinguish numerically from real roots. In ad- 
dition, deeper problems may further arise since the CEP 
is related to a pole with a non-integer critical exponent 
which is generally difficult to describe properly within a 
Pade approximation. 
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FIG. 6. Root distribution in the complex /x/T-plane of a 16 th order Taylor series of the pressure. The corresponding temper- 
atures in the three panels are chosen such that the highest Taylor coefficient cie changes sign while all the remaining Taylor 
coefficients are still positive. 



[1] J. B. Kogut and M. A. Stephanov, Camb. Monogr. Part. 

Phys. Nucl. Phys. Cosmol. 21, 1 (2004). [18 
[2] K. Fukushima and T. Hatsuda, Rept. Prog. Phys. 74, 

014001 (2011). [19 
[3] K. Rajagopal and F. Wilczek (2000), arXiv:hep- 

ph/0011333. 

[4] M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 

(2004) . [20 
[5] C. DeTar and U. M. Heller, Eur. Phys. J. A41, 405 [21 

(2009). [22 

[6] M. Cheng et al, Phys. Rev. D81, 054510 (2010); 
A. Bazavov et al, ibid. D80, 014504 (2009); M. Cheng 
et al, ibid. D74, 054507 (2006). 

[7] S. Borsanyi et al. (Wuppertal-Budapest Collaboration), [23 
JHEP 09, 073 (2010); Y. Aoki et al, ibid. 06, 088 (2009); [24 
Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, Phys. [25 
Lett. B643, 46 (2006). [26 

[8] O. Philipsen, PoS LAT2005, 016 (2006). 

[9] C. Schmidt, PoS LAT2006, 021 (2006). [27 
[10] P. de Forcrand, PoS LAT2009, 010 (2009). 
[11] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, [28 
F. Karsch, E. Laermann, and C. Schmidt, Phys. Rev. 
D66, 074507 (2002); D68, 014507 (2003); C. R. All- [29 
ton, M. Doering, S. Ejiri, S. J. Hands, O. Kaczmarek, 
F. Karsch, E. Laermann, and K. Redlich, D71, 054508 [30 

(2005) . 

[12] R. V. Gavai and S. Gupta, Phys. Rev. D68, 034506 [31 

(2003); D71, 114014 (2005); D78, 114503 (2008). 
[13] C. W. Bernard et al. (MILC Collaboration), Phys.Rev. [32 

D55, 6861 (1997). 
[14] F. Karsch, E. Laermann, and A. Peikert, Phys. Lett. [33 

B478, 447 (2000). 
[15] A. Ali Khan et al. (CP-PACS), Phys. Rev. D64, 074510 [34 

(2001). 

[16] F. Karsch, K. Redlich, and A. Tawfik, Eur. Phys. Jour. [35 

C29, 549 (2003); Phys. Lett. B571, 67 (2003). 
[17] C. Schmidt, PoS CPOD2006, 002 (2006); C. Miao and [36 



C. Schmidt (RBC-Bielefeld) , LAT2008, 172 (2008). 
M. Wagner, A. Walther, and B.-J. Schaefer, Comp. Phys. 
Commun. 181, 756 (2010). 

B. -J. Schaefer, M. Wagner, and J. Wambach, PoS 
CPOD2009, 017 (2009); J. Wambach, B.-J. Schaefer, 
and M. Wagner, Acta Phys. Pol. B Proc. Suppl. 3, 691 
(2010). 

C. Schmidt, PoS PoS CPOD2009, 024 (2009). 
C. DeTar et al, Phys. Rev. D81, 114504 (2010). 

G. A. Baker, Adv. Theo. Phys. 1, 1 (1965); G. A. Baker 
and P. Graves-Morris, Pade approximants; 2nd ed., En- 
cyclopaedia of mathematics and its applications (Cam- 
bridge Univ. Press, Cambridge, 1996). 
M. A. Stephanov, Phys. Rev. D73, 094508 (2006). 
C.-N. Yang and T. Lee, Phys. Rev. 87, 404 (1952). 
F. Karsch, Prog.Part.Nucl.Phys. 62, 503 (2009). 
B.-J. Schaefer, M. Wagner, and J. Wambach, Phys. Rev. 
D81, 074013 (2010). 

B.-J. Schaefer, J. M. Pawlowski, and J. Wambach, Phys. 
Rev. D76, 074023 (2007). 

B.-J. Schaefer and M. Wagner, Progress in Particle and 
Nuclear Physics 62, 381 (2008). 

B.-J. Schaefer and M. Wagner, Phys. Rev. D79, 014018 
(2009). 

V. Skokov, B. Friman, E. Nakano, K. Redlich, and B.-J. 
Schaefer, Phys. Rev. D82, 034029 (2010). 
S. Rofiner, C. Ratti, and W. Weise, Phys. Rev. D75, 
034007 (2007). 

T. K. Herbst, J. M. Pawlowski, and B.-J. Schaefer, 
Phys.Lett. B696, 58 (2011). 

B.-J. Schaefer and J. Wambach, Phys. Rev. D75, 085015 
(2007). 

E. Nakano, B.-J. Schaefer, B. Stokic, B. Friman, and 
K. Redlich, Phys. Lett. B682, 401 (2010). 

F. Karsch, E. Laermann, and C. Schmidt, Phys. Lett. 
B520, 41 (2001). 

Y. Hatta and T. Ikeda, Phys. Rev. D67, 014028 (2003). 



